Marine Aquaculture Suitability Analysis for Oyster and Mussel

Github Repository: https://github.com/jwonyk/marine-aquaculture-selection

Introduction

Exclusive Economic Zones (EEZs) are marine areas extending up to 200 nautical miles from the coast, where the U.S. has special rights to use and manage resources. The growth of marine aquaculture has increased the need for spatial analyses that help identify where different species are most likely to thrive. For shellfish farmers, selecting suitable offshore locations is crucial, which entails determining the proper environmental conditions, such as temperature and depth, that strongly influence survival, growth, overall yield potential, and understanding how to support more sustainable development. This analysis addresses the question:

Which U.S. West Coast EEZ regions provide the most suitable environmental conditions for oyster and mussel aquaculture based on sea surface temperature (SST) and depth?

Data Description

This analysis integrates multiple spatial datasets to evaluate the suitability of aquaculture along the U.S. West Coast (NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1). Sea surface temperature (SST) rasters represent average annual temperatures from 2008–2012 and capture long-term thermal conditions relevant to shellfish growth. These rasters are converted from Kelvin to Celsius to align with biological thresholds reported in the literature.

Bathymetry data (General Bathymetric Chart of the Oceans (GEBCO)) are used to represent ocean depth, which is a key constraint for aquaculture infrastructure and species survival. Depth rasters are resampled and aligned with the SST grid to allow direct cell-by-cell comparison.

Exclusive Economic Zones (EEZs) boundaries are obtained from MarineRegions.org and define the marine areas over which the United States has jurisdiction for resource use up to 200 nautical miles offshore. These boundaries allow suitability to be summarized by management region rather than individual grid cells.

Species-specific temperature and depth thresholds are applied for oysters and mussels to reflect differences in environmental tolerance and habitat requirements are retrieved from SeaLifeBase.

Analysis Setup

We are going to clear the environment and load necessary packages.

# Clear environment and load required packages
rm(list = ls())
pacman::p_load(tidyverse, sf, here, 
               terra, stars, tmap, units,
               ggplot2, viridisLite, dplyr, testthat,
               janitor, knitr, kableExtra, paletteer)

Explore the Data

In this step, I prepare the environmental datasets used to evaluate aquaculture suitability along the U.S. West Coast. I begin by loading multiple years of sea surface temperature rasters, check coordinate reference systems (CRS) to see spatial layers align correctly for later operations, computing their mean to represent long-term thermal conditions, and then converting values from Kelvin to Celsius for ecological interpretation. Depth data are read in and resampled to match the sea surface temperature raster’s resolution, extent, and coordinate reference system, ensuring cell-by-cell comparisons. Both sea surface temperature and depth layers are then reclassified into binary suitable and unsuitable categories based on species-specific environmental thresholds. The resulting suitability rasters are cropped and masked to EEZ boundaries.

# Assign file into lists for average annual sea surface temperature (08-12)
sst_files <- list.files(here("posts", "marine-aquaculture-selection", "data"), 
                      pattern = "average_annual_sst_", 
                      full.names = TRUE)

# Check to make sure all SST files are loading
if(length(sst_files) != 5) {
  
  stop("Fatal error. Please check the file names or path.")
  
} 

# Read in raster objects and shapefile
eez <- st_read(here("posts", "marine-aquaculture-selection", "data", 
                    "wc_regions_clean", "wc_regions_clean.shp"), 
               quiet = TRUE)
depth_rast <- rast(here("posts", "marine-aquaculture-selection", 
                        "data", "depth.tif"))
sst_rast <- rast(sst_files)



# Align EEZ shapefile CRS to SST CRS
if(!identical(st_crs(eez), st_crs(sst_rast))) {
  
  message("EEZ CRS does not match SST CRS! Transforming EEZ to match SST.")
  eez <- st_transform(eez, st_crs(sst_rast))
  
} else {
  
  message("EEZ CRS matches SST CRS.")
  
}

# Match depth raster to SST CRS using `project()`
if(!crs(depth_rast) == crs(sst_rast)) {
  
  message("Depth CRS does not match SST CRS! Projecting depth raster.")
  depth_rast <- project(depth_rast, sst_rast, method = "near")
  
} else {
  
  message("Depth CRS matches SST CRS.")
  
}

Compute Mean SST and Convert to Celsius

Because SST data span multiple years, we average the rasters to obtain a mean temperature surface, which is then converted from Kelvin to Celsius.

# Get the average of SST in Kelvin ()
sst_mean_kelvin <- mean(sst_rast, na.rm = TRUE)

# Convert the average of SST in Celsius
sst_mean_celsius <- sst_mean_kelvin - 273.15

# Rename the column to make sure no confusion later
names(sst_mean_celsius) <- "sst_mean_celsius"

Align Depth Raster to SST Raster

The depth raster must match the SST raster’s resolution and extent to ensure accurate cell-by-cell comparison during suitability calculations.

# Align depth raster to SST
depth_crop <- crop(depth_rast, sst_mean_celsius)

# Resample to match SST resolution (use the nearest neighbor)
depth_resample <- resample(depth_crop, 
                           sst_mean_celsius, 
                           method = "near")

# Rename the column to make sure no confusion later
names(depth_resample) <- "depth_resample"

Reclassifying & Combine for Suitability

We convert continuous SST values into a binary suitability layer based on oyster temperature thresholds from the literature. Similarly, depth values are reclassified to identify areas within the allowable depth range for oyster growth. Multiplying the two binary layers identifies grid cells that meet both environmental requirements.

# Oyster thresholds
oyster_sst_min <- 11
oyster_sst_max <- 30
oyster_depth_min <- -70
oyster_depth_max <- 0

# Reclassify SST suitability
# 0 = unsuitable & 1 = suitable
oyster_sst_suitable <- lapp(sst_mean_celsius, 
                            fun = function(x)
                              ifelse(x >= 11 & x <= 30, 1, 0))


# Reclassify depth
oyster_depth_suitable <- lapp(depth_resample,
                              fun = function(x)
                                ifelse(x >= -70 & x <= 0, 1, 0))

# Combine classified SST and depth by multiplying them both
oyster_cells_suitable <- oyster_sst_suitable * oyster_depth_suitable

Crop and Mask Suitability to EEZ Boundaries

We restrict the area of interest to U.S. West Coast Exclusive Economic Zones (EEZs), ensuring suitability is evaluated only within jurisdictional waters.

# Convert EEZ object to SpatVector to use as `terra` ojbect
eez_vector <- vect(eez)

# Crop suitability raster to EEZ bounding box (reduce processing)
oyster_crop <- crop(oyster_cells_suitable, eez_vector)

# Mask cells outside the EEZ
oyster_eez <- mask(oyster_crop, eez_vector)

Compute Suitable Cells to Area

Each raster cell’s area is computed and summed within each EEZ to quantify total suitable habitat in square kilometers.

# Compute per-cell area (in km^2)
cell_area_km2 <- cellSize(oyster_eez, unit = "km")

# Multiply area raster by suitable, 0 vs 1, only return suitable cells
oyster_area <- oyster_eez * cell_area_km2

# Rename the column to make sure no confusion later
names(oyster_area) <- "suitable_area_km2"

Join Suitability Results to EEZ Polygons

We join calculated suitable area values back to EEZ polygons so they can be visualized on a choropleth map.

# Rename ID column to easier access
eez <- eez %>% 
  rename(eez_id = rgn_id)

# Convert into SpatVector after name update
eez_vector <- vect(eez)

# Rasterize EEZ polygons to match raster grid before `zonal()`
eez_id_raster <- rasterize(eez_vector, oyster_area, field = "eez_id")

Calculate and Summarize Suitable Area by EEZ

We sum all suitable raster cells within each EEZ and join these results to the EEZ polygons.

# Add up suitable area within each EEZ zone
oyster_area_eez <- zonal(oyster_area, eez_id_raster,
                         fun = "sum",
                         na.rm = TRUE) %>% 
  as.data.frame()

# Join with EEZ polygons for mapping purposes
eez_oyster_join <- eez %>% 
  left_join(oyster_area_eez, by = "eez_id")

# Create table of suitable area ranked by EEZ
eez_oyster_join %>% 
  st_drop_geometry() %>% 
  arrange(desc(suitable_area_km2)) %>% 
  kable("html", 
        col.names = c("Region", "Region Code", "Area(m^2)",
                      "EEZ ID", "EEZ Area (km^2)", "Suitable Area (km^2)"),
        caption = "Suitable Oyster Habitat Area by EEZ (km^2)",
        align = "c") %>% 
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover"))

Suitable Oyster Habitat Area by EEZ (km^2)

Region	Region Code	Area(m^2)	EEZ ID	EEZ Area (km^2)	Suitable Area (km^2)
Central California	CA-C	202738329147	3	202738.33	4940.0422
Southern California	CA-S	206860777840	4	206860.78	4221.3919
Washington	WA	66898309678	5	66898.31	3313.1592
Oregon	OR	179994061293	1	179994.06	1578.9688
Northern California	CA-N	164378809215	2	164378.81	454.2957

Map Oyster Aquaculture Suitability

Oysters require relatively warm sea surface temperatures and moderate depths, making their potential habitat strongly influenced by latitudinal temperature gradients along the West Coast. The map below visualizes how these environmental constraints translate into regional differences in oyster aquaculture suitability across EEZs, highlighting areas where conditions are most favorable for oyster growth and infrastructure deployment.

# Extract paleteer colors
teal_palette <- as.character(paletteer_d("ggsci::teal_material"))
eez_map <- eez_oyster_join %>% 
  filter(!is.na(suitable_area_km2)) 

# Map suitable area by EEZ
tm_tiles("Esri.OceanBasemap") +
tm_shape(eez_map) +
  tm_polygons(
    col = "suitable_area_km2",
    palette = colorRampPalette(teal_palette)(5),
    title = "Suitable Area (km^2)") +
  tm_layout(main.title = "Oyster Aquaculture Suitability by EEZ",
            main.title.size = 1,
            frame = FALSE,
            legend.outside = TRUE,
            legend.outside.position = "right",
            outer.margins = c(0, 0, 0, 0),
            inner.margins = c(0, 0, 0, 0.02)) +
  tm_compass(position = c("right", "top")) +
  tm_scalebar(breaks = c(0, 100, 200, 300),
              position = c("left", "bottom"))

Reflection - Oyster Aquaculture Suitability

The oyster aquaculture suitability analysis demonstrates how combining raster-based environmental data with vector-based EEZ boundaries reveals spatial patterns relevant to marine resource planning. By averaging five years of SST data and resampling bathymetry, this analysis identifies regions where environmental conditions align with oyster growth requirements. Warmer southern EEZs exhibit substantially larger suitable areas than northern regions, highlighting how temperature and depth jointly constrain oyster aquaculture potential along the West Coast.

Apply Workflow Function

To show generalizability, the suitability workflow is applied to another species with different environmental thresholds.

# General Function to calculate aquaculture suitability

compute_suitability <- function(sst_min, sst_max,
                                depth_min, depth_max,
                                species_name) {
  
# 1. Reclassify SST
sst_suitable <- lapp(sst_mean_celsius,
                     fun = function(x) 
                       ifelse(x >= sst_min & x <= sst_max, 1, 0))

# 2. Reclassify Depth
depth_suitable <- lapp(depth_resample,
                       fun = function(x) 
                         ifelse(x >= depth_min & x <= depth_max, 1, 0))

# 3. Combine Suitability
suitable_cells <- sst_suitable * depth_suitable

# 4. Crop and Mask to EEZ
suitable_crop <- crop(suitable_cells, eez_vector)
suitable_eez <- mask(suitable_crop, eez_vector)

# 5. Convert Suitability to Area in km^2
cell_area <- cellSize(suitable_eez, unit = "km")
suitable_area <- suitable_eez * cell_area
names(suitable_area) <- "suitable_area_km2"

# 6. Rasterize EEZ and Compute Zonal Area
eez_id_raster <- rasterize(eez_vector, suitable_area, field = "eez_id")

zone_df <- zonal(suitable_area, eez_id_raster,
                 fun = "sum", na.rm = TRUE) %>% 
  as.data.frame()

# 7. Join results
eez_join <- eez %>% 
  left_join(zone_df, by = "eez_id")

# 8. Map
teal_palette <- as.character(paletteer_d("ggsci::teal_material"))

map <- tm_tiles("Esri.OceanBasemap") +
  tm_shape(eez_join) +
    tm_polygons(col = "suitable_area_km2",
                palette = colorRampPalette(teal_palette)(20),
                title = "Suitable Area (km²)") +
    tm_layout(main.title = paste("Aquaculture Suitability for",
                                 species_name),
              main.title.size = 1.0,
              legend.outside = TRUE,
              frame = FALSE) +
    tm_compass(position = c("right", "top")) +
    tm_scalebar(breaks = c(0, 100, 200, 300),
                position = c("left", "bottom"))
return(map)}

Map Mussel Aquaculture Suitability

Mussels are adapted to cooler waters and shallower depths, leading to a spatial suitability pattern that differs markedly from that of oysters. The map below illustrates how these contrasting environmental preferences shift suitable habitat northward, revealing regions where mussel aquaculture may be more viable under existing temperature and depth conditions.

# SST 9 to 18 degrees Celsius, depth 0 to 20 meters
mussle_map <- compute_suitability(
  sst_min = 9,
  sst_max = 18,
  depth_min = -20,
  depth_max = 0,
  species_name = "Mussels (A. californiensis)")

mussle_map

Reflection - Mussel Aquaculture Suitability

The mussel suitability assessment shows different spatial patterns compare to oysters, emphasizing how species-specific environmental requirements could lead to different aquaculture potential across EEZs. Mussels favor cooler temperature and shallow depths and shows that more suitable northward than southern regions - sustainability is highest in Northern California and Washington.

Conclusion

This analysis demonstrates how we use spatial data science to evaluate the sustainability of aquaculture across large marine regions. By combining sea surface temperature data, bathymetry, raster reclassification, and Zonal statistics, we can identify locations that hit all the requirements for the environmental conditions specific to the aquaculture species. The results show that oysters prefer the warm southern waters, whereas mussels are better suited to the cooler northern regions of EEZs. The findings could improve productivity and reduce ecological risks.

Future Directions

While this analysis highlights broad spatial patterns in aquaculture suitability, several limitations should be addressed in future work. Additional environmental variables, such as salinity, chlorophyll concentration, ocean currents, and wave exposure, play important roles in shellfish growth but are not included here. Socioeconomic and regulatory constraints also influence where aquaculture development is feasible, as not all environmentally suitable locations are accessible or permitted. Incorporating higher-resolution datasets and seasonal variability could further refine suitability estimates and support more detailed marine spatial planning.

Reference

GEBCO Compilation Group. (2022). GEBCO_2022 Grid. Retrieved from https://www.gebco.net

Gentry, R. R., Froehlich, H. E., Grimm, D., Kareiva, P., Parke, M., Rust, M., Gaines, S. D., & Halpern, B. S. (2017). Mapping the global potential for marine aquaculture. Nature Ecology & Evolution, 1, 1317–1324.

Hall, S. J., Delaporte, A., Phillips, M. J., Beveridge, M., & O’Keefe, M. (2011). Blue Frontiers: Managing the Environmental Costs of Aquaculture. The WorldFish Center.

Marine Regions. (2025). Exclusive Economic Zones data. Retrieved from https://www.marineregions.org

National Oceanic and Atmospheric Administration. (2008–2012). 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1. Retrieved from https://www.ncei.noaa.gov

SeaLifeBase. (2025). Adula californiensis species summary. Retrieved from https://www.sealifebase.ca/summary/Adula-californiensis.html